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ABSTRACT 

Convex optimization problems are common in hyperspectral unmix- 
ing. Examples are the constrained least squares (CLS) problem used 
to compute the fractional abundances in a linear mixture of known 
spectra, the constrained basis pursuit (CBP) to find sparse {i.e., 
with a small number of terms) linear mixtures of spectra, selected 
from large libraries, and the constrained basis pursuit denoising 
(CBPDN), which is a generalization of BP to admit modeling er- 
rors. In this paper, we introduce two new algorithms to efficiently 
solve these optimization problems, based on the alternating direction 
method of multipliers, a method from the augmented Lagrangian 
family. The algorithms are termed SUnSAL {sparse unmixing by 
variable splitting and augmented Lagrangian) and C-SUnS AL {con- 
strained SUnSAL). C-SUnSAL solves the CBP and CBPDN prob- 
lems, while SUnSAL solves CLS as well as a more general version 
thereof, called constrained sparse regression (CSR). C-SUnSAL 
and SUnSAL are shown to outperform off-the-shelf methods in 
terms of speed and accuracy. 

1. INTRODUCTION 

Hyperspectral unmixing (HU) is a source separation problem with 
applications in remote sensing, analytical chemistry, and other areas 
ll0llll|[T2l . Given a set of observed mixed hyperspectral vectors, 
HU aims at estimating the number of reference spectra (the end- 
members), their spectral signatures, and their fractional abundances, 
usually under the assumption that the mixing is linear 1 1011 121 . 

Unlike in a canonical source separation problem, the sources in 
HU {i.e., the fractional abundances of the spectra/materials present 
in the data) exhibit statistical dependency 1151 . This characteristic, 
together with the high dimensionality of the data, places HU beyond 
the reach of most standard source separation algorithms, thus foster- 
ing active research in the field. 

Most HU methods can be classified as statistical or geometri- 
cal. In the (statistical) Bayesian framework, all inference relies on 
the posterior probability density of the unknowns, given the observa- 
tions. According to Bayes' law, the posterior probability density re- 
sults from two factors: the observation model (the likelihood), which 
formalizes the assumed data generation model, possibly including 
random perturbations such as additive noise; the prior, which may 
impose natural constraints on the endmembers {e.g., nonnegativity) 
and on the fractional abundances {e.g., belonging to the probability 
simplex, since they are relative abundances), as well as model spec- 
tral variability (5lfT4l[T6). 

Geometrical approaches exploit the fact that, under the linear 
mixing model, the observed hyperspectral vectors belong to a sim- 
plex set whose vertices correspond to the endmembers. Therefore, 
finding the endmembers amounts to identifying the vertices of that 
simplex fT8][T6ll3l[T3ll20l[Tl. 



Sparse regression is another direction recently explored for HU 
]9j, which has connections with both the statistical and the geomet- 
rical frameworks. In this approach, the problem is formulated as 
that of fitting the observed (mixed) hyperspectral vectors with sparse 
{i.e., containing a small number of terms) linear mixtures of spec- 
tral signatures from a large dictionary available a priori. Estimating 
the endmembers is thus not necessary in this type of methods. No- 
tice that the sparse regression problems in this context are not stan- 
dard, as the unknown coefficients (the fractional abundances) sum to 
one (the so-called abundance sum constraint - ASC) and are non- 
negative {abundance non-negativity constraint - ANC). These prob- 
lems are thus referred to as constrained sparse regression (CSR). 

Several variants of the CSR problem can be used for HU; some 
examples follow. In the classical constrained least squares (CLS), 
the fractional abundances in a linear mixture of known spectra are 
estimated by minimizing the total squared error, under the ANC and 
the ASC. Although no sparseness is explicitly encouraged in CLS, 
under some conditions (namely positivity of the spectra) it can be 
shown that the solutions are indeed sparse f2). Constrained basis 
pursuit (CBP) is a variant of the well-known basis pursuit (BP) cri- 
terion 13) under the ANC and the ASC; as in BP, CBP uses the i\ 
norm to explicitly encourage sparseness of the fractional abundance 
vectors. Finally, constrained basis pursuit denoising (CBPDN) is a 
generalization of CBP that admits modeling errors (e.g., observation 
noise). 

1.1. Contribution 

In this paper, we introduce a class of alternating direction algorithms 
to solve several CSR problems (namely CLS, CBP, and CBPDN). 
The proposed algorithms are based on the alternating direction 
method of multipliers (ADMM) |8 7 6|, which decomposes a diffi- 
cult problem into a sequence of simpler ones. Since ADMM can be 
derived as a variable splitting procedure followed by the adoption of 
an augmented Lagrangian method to solve the resulting constrained 
problem, we term our algorithms as SUnSAL (spectral unmixing by 
splitting and augmented Lagrangian) and C-SUnSAL (constrained 
SUnSAL). 

The paper is organized as follows. Section|2]introduces notation 
and formulates the optimization problems. Section 3 reviews the 
ADMM and the associated convergence theorem. Section 4 intro- 
duces the SUnSAL and C-SUnSAL algorithms. Section 5 presents 
experimental results, and Section 6 ends the paper by presenting a 
few concluding remarks. 

2. PROBLEM FORMULATION: CSR, CLS, CBP, CBPDN 

Let A € R fcxn denote a matrix containing the n spectral signatures 
of the endmembers, x £ R n denote the (unknown) fractional abun- 
dance vector, and y £ R fc be an (observed) mixed spectral vector. 



In this paper, we assume that A is known; this is the case the CSR 
approach [ 9 1, where it is a library with a large number of spectral sig- 
natures, thus usually n > k. Matrix A can also be the output of an 
endmember extraction algorithm, in which case usually n <C k. The 
key advantage of the CSR approach is that it avoids the estimation 
of endmembers, quite often a very hard problem. 
The general CSR problem is defined as 

L(x) 

, *• , 

(P CSR ): min (l/2)||Ax-y||! + A||x||i (1) 

X 

subject to: x > 0, l T x = 1, 

where ||x||2 and ||x||i denote the I2 and t\ norms of x, respectively, 
A > is a parameter controlling the relative weight between the £2 
and £x terms, 1 denotes a column vector of 1 's, and the inequality 
x > is to be understood in the componentwise sense. The con- 
straints x > and l T x = 1 correspond to the ANC and ASC, re- 
spectively. The CLS problem corresponds to Pcsr with A = 0. The 
CBP and CBPDN problems are also equivalent to particular cases of 
Pcsr, as stated next. 

The CBP optimization problem is 

(P CB p): min ||x||i (2) 
subject to: Ax = y, x > 0, l T x = 1, 

Notice that Pcbp corresponds to Pcsr with A — > do. 
The CBPDN optimization problem is 

(Pcbpdn): min ||x||i (3) 

X 

subject to: || Ax - y|| 2 < 5, x > 0, l T x = 1. 

Problem Pcsr ls equivalent to Pcbpdn m the sense that for any 
choice of S for which Pcbpdn is feasible, there is a choice of X for 
which the solutions of the two problems coincide j 1 911 . Finally, no- 
tice that Pcbp corresponds to Pcbpdn with 5 = 0. 

3. THE ADMM 

Consider an unconstrained problem of the form 

min / 1 (x) + / 2 (Gx), (4) 

where fx : R™ ft, f 2 : R p ft, and G G R px ". The ADMM 
(U|7][8], me tool in this paper, is as shown in Fig.[T] The follow- 
ing is a simplified version of a theorem of Eckstein and Bertsekas 
stating convergence of ADMM. 

Theorem 1 (| 6 1) Lef G have full column rank and fx , fi be closed, 
proper, and convex. Consider arbitrary /j, > and uu,do G W. 
Consider three sequences {xfc G R™, k — 0,1,...}, {life G 
R p , k = 0, 1, ...}, and{d fe G R p , k = 0, 1, ...} that satisfy 

x fc+ i = argmin/i(x) + — IIGx— u*— cMij (5) 
x 2 

u k+1 = argmin/ 2 (u) + -||Gx fc+ i-u-difc||2 (6) 

u 2 

d k+1 = d fe - (Gx fc+ i - Ufc+i). (7) 

Then, if @) has a solution, the sequence {xfe} converges to it; other- 
wise, at least one of the sequences {u k } or {d k } diverges. 



Algorithm ADMM 

1 . Set k = 0, choose /i > 0, uo , and do . 

2. repeat 

3. x fe+1 G argmin x /i(x) + f ||Gx - u fc - dfc ||| 

4. u fc+1 G argmm u / 2 (u) + ^||Gx fc+1 - u - d fc ||| 

5. d fe+ i 4- d fe - (Gx fe+ i - Ufc+i) 

6. k 4-k + l 

7. until stopping criterion is satisfied. 

Fig. 1. The alternating direction method of multipliers (ADMM). 

4. APPLICATION OF ADMM 

In this section, we specialize the ADMM to each of the optimization 
problems stated in Section[2] 

4.1. ADMM for CSR and CLS: the SUnSAL Algorithm 

We start by writing the optimization Pcsr in the equivalent form 

min (1/2) || Ax - y||l + A||x||i + (l T x) + t K » (x), (8) 

where ls is the indicator function of the set S (i.e., ts(x) = if 
x G S and ts(x) = cxj if x ^ S). We now apply the ADMM using 
the following translation table: 

/ 1 (x) = i||Ax-y||i + t{1} (l T x) (9) 

/ 2 (x)sA||x||i + i. R «(x) (10) 

GeI. (11) 

With the current setting, step 3 (Fig. 1) of the ADMM requires 
solving a quadratic problem with linear equality constraint, the solu- 
tion of which is 

Xfc+i <- B _1 w - C(l T B _1 w - 1) (12) 

where 

B = A T A + ill (13) 

C ee B _1 1(1 T B _1 1) _1 (14) 

w ee A T y + L i(u k +d fc ). (15) 

Step 4 of the ADMM (Fig. 1) is simply 

Ufe+i <- argmin (l/2)||u - u k \\l + (X/ fi)\\u\\x + t K «(u) (16) 

where v k ee Xfe + i — dfe. Without the term ign , the solution of l |16t 
would be the well-known soft threshold (4): 

Ufc+i 4r- S0ft(l/fc, (17) 

A straightforward reasoning leads to the conclusion that the effect of 
the ANC term is to project onto the first orthant, thus 

Ufe + i <- max{0, soft(j/fe, X/li)}, (18) 

where the maximum is to be understood in the componentwise sense. 

Fig. [2] shows the SUnSAL algorithm, which is obtained by re- 
placing lines 3 and 4 of ADMM by l !12t and dl8| l. respectively. 

The objective function {8]l is proper, convex, lower semi- 
continuous, and coercive, thus it has a non-empty set of minimizers 



Algorithm SUnSAL 

1. Set k = 0, choose fi > 0, Urj, and do 
repeat 

w +- A T y + /i(u fc + dfc) 
x fc+i 



Algorithm C-SUnSAL 



B^w - C(l T B- 1 w - 1) 



<- *k+l - dfc 
Ufc+i <- max{0,soft(i/fc, A//i)} 
dfc+i <- d fe - (Xfc + i - u fe+1 ) 
k «- fc + 1 
until stopping criterion is satisfied. 



Fig. 2. Spectral unmixing by variable slitting and augmented La- 
grangian (SUnSAL). 

(see 1191 . for definitions of these convex analysis concepts). Func- 
tions /i and / 2 in l[9) and d 1 Ofc are closed and G = I is obviously 
of full column rank, thus Theorem 1 can be invoked to ensure 
convergence of SUnSAL. 

Concerning the computational complexity, we refer that, in hy- 
perspectral applications, the rank of matrix B is no larger that the 
number of bands, often of the order of a few hundred, thus B 1 can 
be easily precomputed. The complexity of the algorithm per itera- 
tion is thus 0(n 2 ), corresponding to the matrix-vector products. 

To solve the CLS problem, we simply run SUnSAL with A = 0. 
Moreover, the ANC (non-negativity constrain) can be turned off by 
using d!7t instead of (18) and the ASC (sum to one constraint) can 
be deactivated by setting C = in H2\ . 

4.2. ADMM for CBP and CBPDN: the C-SUnSAL Algorithm 

Given that the CBP problem corresponds to CBPDN with S = 0, we 
address only the latter. Problem Pcbpdn ls equivalent to 



min ||x||i + i B(yij )(Ax) + i{i}(l T x) + t R+ (x) 



(19) 



where B(y,S) = {z : ||z — y || 2 < S} is a radius-5 closed ball 
around y. To apply the ADMM we use the following definitions: 

/i(x) = l{ i } (l T x) (20) 
/a(u) = i B ( y ,j)(ui) + A||u 2 ||i + t R » (u 2 ) (21) 
G = [A T I] T . (22) 

where u = [ui T u^] T . With the above definitions, the solution of 
line 3 of ADMM (Fig. 1), a quadratic problem with linear equality 
constraints, is 



Xfc+l 



B _1 w-C(l T B~ 1 w- 1) 



where 



B = A A + I 

C = B _1 1(1 T B _1 1) _1 

w = A T (ui,fc + di.fe) + (u 2 ,fc + d 2 ,fc 



(24) 
(25) 
(26) 



Because the variables ui and u 2 are decoupled, line 4 of 
ADMM (Fig. 1) consists in solving two separate problems, 

ui,hi £ argmin (l/2)||u-i/i lfe ||2 + ts(y,a)(u) (27) 

u 2 ,fc+i £ argmin (1/2) ||u - i/2,fc||| + (A/^)||u||i + m« (u) 

(28) 



1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 

12. 



Set k <— 0, choose > 0, Ui n> dio, u 2 .o, and d 2j o- 
repeat 

w <- A T (u ljfc + d 1>fe ) + (u 2 ,/b + d a ,fc) 
x fe+ i <- B^w - C(l T B- 1 w - 1) 
"l.fc <- Ax fc+1 - d l fe 

ui, fc+ i <- ^(y.^C^i.fe) 

"2,k <— Xfe + i — d 2t . 
u 2,fc+i <- max{0, soft(i/ 2 ,fc, Vm)} 
di,fc+i <- d life - (Ax fc+1 - u ljfe+1 ) 
d 2 ,fc+i <— d 2fe — (x fc+1 — u 2:ft+1 ) 
k <- fc + 1 
until stopping criterion is satisfied. 



Fig. 3. Constrained spectal unmixing by variable slitting and aug 
mented Lagrangian (C-SUnSAL). 



where 

i/i,* = Ax fe+ i - d ljfc (29) 
fa, fc = Xfc+i - d 2 ,fe. (30) 

The solution of l !27t is the projection onto the ball B(y, S), given by 



"ik, \Wi,k - y II 2 < s 



ll"i,k-y||a 

Similarly to j 1 8b . the solution of d28t is given by 

U2,fc+i <- max{0, soft(i/ 2 ,fe, A//x)}. 



(31) 



(32) 



Fig. [3] shows the C-SUnSAL algorithm for CBPDN, which re- 
sults from replacing line 3 of ADMM (Fig. 1) by J23l > and line 4 of 
ADMM by (alHSD. As mentioned above, C-SUnSAL can be used 
to solve the CBP problem simply by setting S — 0. As in SUnSAL, 
the ANC and/or the ASC can be deactivated trivially. 

The objective function l |19l > is proper, convex, lower semi- 
continuous, and coercive, thus it has a non-empty set of minimizers. 
Functions /1 and / 2 in l |20| > and ( 1211 ) are closed and G in d22t is 
obviously of full column rank, thus Theorem 1 can be invoked to 
ensure convergence of C-SUnSAL. 

Concerning the computational complexity, the scenario is simi- 
lar to that of SUnSAL, thus complexity of C-SUnSAL is 0(n 2 ) per 
iteration. 

5. EXPERIMENTS 



(23) We now report experimental results obtained with simulated data 
generated according to y = Ax + n, where n £ R fc models addi- 



tive perturbations. In hyperspectral applications, these perturbations 
are mostly model errors dominated by low-pass components. For 
this reason, we generate the noise by low-pass filtering samples of 
zero-mean i.i.d. Gaussian sequences of random variables. We define 
the signal-to-noise ratio (SNR) as 



SNR= 10 log! 



E[||Ax|| 

E[||n|| = 



The expectations in the above definition are approximated with sam- 
ple means over 10 runs. The original fractional abundance vectors 
x are generated in the following way: given s, the number of non- 
zero components in x, we generate random samples uniformly in the 
(s — 1)— simplex and distribute randomly these s values among the 



Table 1. RSNR values and execution times for the Gaussian library 
defined in the text (average over 10 runs). 
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Table 2. RSNR values and execution times for the USGS library 
(average over 10 runs). 
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components of x. We considered two libraries (i.e., matrices A): 
a 200 x 400 matrix with zero-mean unit variance i.i.d. Gaussian 
entries and a 224 x 498 matrix with a selection of 498 materials 
(different mineral types) from the USGS library denoted splibOffl. 

As far as we know, there are no special purpose algorithms for 
solving the CSR, CBP, and CBPDN problems. Of course these 
are canonical convex problems, thus they can be tackled with stan- 
dard convex optimization techniques. Namely, the CLS, which is a 
particular case of CSR, can be solved with the MATLAB function 
lsqnonneg, which we use as baseline in our comparisons. 

Tables 1 and 2 report reconstruction SNR (RSNR), defined as 

where x is the estimated fractional abundance vector, and execution 
times, for the two libraries referred above. The lsqnonneg is run 
with its default options. SUnSAL and C-SUnSAL run 200 iterations, 
which was found to be more than enough to achieve convergence. 

We highlight the following conclusions: (a) the proposed algo- 
rithms achieve higher accuracy in about two orders of magnitude 
shorter time. This is a critical issue in imaging application where an 
instance of the problem has to be solved for each pixel; (b) the lower 
accuracy obtained with the USGS matrix is due to the fact that the 
spectral signatures are highly correlated resulting in a much harder 
problem than with the Gaussian matrix. 

6. CONCLUDING REMARKS 

In this paper, we introduced new algorithms to solve a class of opti- 
mization problems arising in spectral unmixing. The proposed algo- 
rithms are based on the alternating direction method of multipliers, 
which decomposes a difficult problem into a sequence of simpler 
ones. We showed that sufficient conditions for convergence are sat- 
isfied. In limited set of experiments, the proposed algorithms were 
shown to clearly outperform an off-the-shelf optimization tool. On- 
going work includes a comprehensive experimental evaluation of the 
proposed algorithms. 



'http://speclab.cr.usgs.gov/spectral.lib06 
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